Striped spin liquid crystal ground state instability of kagome antiferromagnets 
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The Dirac spin liquid ground state of the spin 1/2 Heisenberg kagome antiferromagnet has po- 
tential instabilities [HI]- This has been suggested as the reason why it is not strongly supported 
in large-scale numerical calculations 5.. However, previous attempts to observe these instabilities 
have failed. We report on the discovery of a projected BCS state with lower energy than the pro- 
jected Dirac spin liquid state which provides new insight into the stability of the ground state of 
the kagome antiferromagnet. The new state has three remarkable features. First, it breaks both 
spatial symmetry in an unusual way that may leave spinons deconhned along one direction. Second, 
it breaks the U(l) gauge symmetry down to Z2- Third, it has the spatial symmetry of a previously 
proposed "monopole" suggesting that it is an instability of the Dirac spin liquid. The state described 
herein also shares a remarkable similarity to the distortion of the kagome lattice observed at low Zn 
concentrations in Zn-Paratacamite suggesting it may already be realized in these materials. 



Many potential ground states have been suggested 
for the spin 1/2 kagome Heisenberg antiferromagnet; 
these include magnetic ordering in one of several spin 
patterns [B], a valence bond crystal with a 36 site unit 
cell[7] or a 12 site unit cell[S], a chiral spin liquid[7J [5], 
several kinds of Z 2 spin liquids EH and an algebraic 
spin liquid [U |3]. Amongst algebraic spin liquids, it was 
shown [3 that the [/(l)-Dirac state had lowest variational 
energy. This latter state is characterized by its graphene- 
like Dirac fcrmionic spinons interacting with a Maxwell- 
like photon that characterizes singlet excitations. Many 
of the other suggestions for the kagome antiferromag- 
net ground state can be characterized as instabilities of 
this Dirac spin liquid; these instabilities have been cata- 
logued in ref. |31 However, despite the many instabilities, 
all variational studies we are aware of that include the 
[/(l)-Dirac state for some choice of variational parame- 
ters have found it remarkably stable [2ti4) Hi] fl2] . 

Recently, large scale DMRG calculations [SJ [T31 [H] 
have produced strong numerical evidence that the true 
ground state of the kagome antiferromagnet is a Z 2 spin 
liquid. In particular, they have found the lowest known 
variational state with energies comparable to exact di- 
agonalization on small clusters, a gap to all excitations 
and topological entanglement entropy expected of a Z 2 
spin liquid. These results are quite surprising in light of 
prior studies of the U(l)-Dirac state and suggest that it 
is unstable. 

Perhaps the simplest explanation of the DMRG results, 
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FIG. 1. The n.n. s-s correlations for the projected Dirac Spin 
liquid state with twisted boundary conditions (left) and the 
lowest energy projected mean field state (right). The color 
indicates sign and the linewidth measures the magnitude. 



as proposed in Ref. [TUJ is that the Dirac-like spinons 
pair up and form a superfluid-like Z 2 spin liquid. With 
this goal in mind, these authors catalog 20 such possi- 
ble states and give explicit instructions for constructing 
14 of them. Unfortunately, despite searches for such a 
state[T2], no Z 2 state has been shown to be energetically 
favorable to the U{1) Dirac state at the level of varia- 
tional projected BCS wave functions. 

In this letter, we revisit the projected BCS (PBCS) 
variational wave function problem on the kagome lattice 
and have discovered an instability of the Dirac state. Our 
approach, building on Refs. [15] and 1161 is to optimize 
over the entire set of time-reversal invariant projected 
BCS wave functions which are parameterized by a sym- 
metric pairing matrix 4>ij. These states encompass the 
U(l) and Z2 spin liquids, general valence bond solids, 
and a number of other instabilities. We find states with 
lower energy than the Dirac state on clusters of up to 192 
sites. On the 48 site cluster in particular, a carefully op- 
timized wave function, as shown in Fig. [T] demonstrates 
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that the new state breaks lattice symmetries, doubling 
the unit cell with a wave vector that connects the two 
Dirac nodes. However it preserves a C% rotational sym- 
metry about several lattice points giving it a distinct one 
dimensional character. It also lies very close to the Dirac 
state with similar fluxes through the hexagons and bow 
ties. For these reasons, we believe it should belong to one 
of the leading instabilities proposed in Ref. [4] and have 
discovered that the most likely candidate is the time- 
reversal symmetric "w-monopole" . We also find that this 
state cannot be obtained by optimizing within the class 
of wave functions that do not admit pairing and therefore 
likely has Z2 global gauge symmetry. For this reason and 
because of its Ci symmetry and one dimensional charac- 
ter, we have dubbed this state a Z 2 striped (or smectic) 
spin liquid crystal. Finally, we discuss a connection be- 
tween our state and Zn-Paratacamite with Zn concentra- 
tions x below x = 1/3 whose observed distorted kagome 
layers[17] has the same symmetry as Fig. [ijb). 

Our objective is to find the ground state of the nearest 
neighbor spin 1/2 kagome antiferromagnet constrained 
to the set of Gutzwiller projected variational BCS wave 
functions 

l*> =£s=i/2|*o) (1) 

where * (-R) = (-R|*o) = det M{R), M tJ (R) = 
4>(rif,fji) is the BCS pairing amplitude and Ps=i/2 
projects these states onto the physical singly occupied 
Hilbert space of spin wave functions. This set of varia- 
tional ansatz includes projected Slater determinants such 
as the U(X) Dirac spin liquid [3 . The variational search 
is performed by minimizing the ground state energy 

E={*\H\9), H = Y,S?S?. (2) 

(ij) 

of the Heisenberg spin Hamiltonian on the kagome lat- 
tice using as parameters all n(n — l)/2 real values of the 
pairing function (f>(fi,fj) (building on Refs. IT51 and ITS)) . 
To avoid only finding local minima in this energy land- 
scape, we choose to start the optimization in many qual- 
itatively different starting points by making use of the 
well developed projective symmetry group classification 
developed in the literature[cO[Il[TOJ[TH]. In particular, we 
will choose mean field parameters following Ref. [TO] and 
derive a pairing matrix from these following supplemen- 
tal materials (SM) S-I and STL This allows us to start 
from 14 different spin liquid wave functions. Fig. [2^a) 
shows prototypical optimization traces of the energies for 
each spin liquid state on the 4x4 lattice, labeled follow- 
ing Table II of Ref. 10. Many states lie below the energy 
of the Dirac spin liquid and the minimal energy state we 
find is at -0.430520 ± 5 x 10~ 6 per site (see SM S-III for 
its pairing matrix) , significantly below the Dirac spin liq- 
uid energy of 0.42938 ±4 x 10 -5 per site showing that the 
Dirac spin liquid is not the most stable projected mean 
field state. 
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FIG. 2. (a): Prototypical stochastic optimization traces start- 
ing from the 14 distinct PSG ansatz defined in Table II of 
Ref. [10] and optimizing the pairing function (f>. Inset shows 
magnified version, (b): Optimization of only single particle 
orbitals starting from perturbed Dirac states where X labels 
the strength of the perturbation. 

This is particularly surprising in light of previous stud- 
ies that find the Dirac spin liquid to be stable against 
many instabilities^, 3, 4][T2][T9]. Such calculations failed 
to find this state because they only considered short- 
range mean field Hamiltonians, a limitation we avoid by 
optimizing a general pairing function that allows both 
long-range hopping terms and spatial symmetry breaking 
in the fermion Hamiltonian. As shown in Fig. [2jb) and 
discussed in more detail below, our approach in conjunc- 
tion with enforcing U(l) symmetry does not lower the 
energy. But, additionally allowing the U(l) gauge sym- 
metry to break down to Z2 does lower the energy (left 
side of Fig. [2]) . Consequently, the gauge symmetry break- 
ing down to Z2 is intricately linked to spatial symmetry 
breaking and also possibly long-range hopping/pairing 
terms in the Hamiltonian. 

Let us then focus on the nature of our lowest energy 
state. Most strikingly we find that the state breaks trans- 
lational symmetry doubling the unit cell. Fig. [I] (right) 
shows the deviation of the spin-spin correlation function 
on nearest neighbors from its average value, in contrast 
with the U(l) Dirac spin liquid shown in Fig. [I] (left). 
The symmetry breaking seen in the Dirac spin liquid is 
a finite size artifact coming from twisted boundary con- 
ditions in one direction. In our optimized state we find 
breaking of the translational symmetry and a doubling 
of the unit cell leading us to conclude that the minimal 
projected mean field state is not an isotropic spin liquid. 

To further verify that there is not a nearby isotropic 
spin liquid state to our optimal state, we stochastically 
modify the pairing function, starting from a low energy 
state, in small increments searching for a true symmetric 
spin liquid. In particular, we make a random stochastic 
change to the pairing function accepting it only if we 
have moved closer to a spin liquid state (defined by the 
deviations of the unprojected (Si ■ Sj)). In this way, we 
search for the "closest" spin liquid. When we run this 
procedure, the spin liquid state it approaches has the 
energy of the Dirac spin liquid (and hence is presumably 
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the Dirac spin liquid). This leads us to believe that their 
is no additional nearby isotropic spin liquid. 

Having found a broken symmetry state, we would like 
to now understand whether it retains characteristics rem- 
iniscent of any spin liquid phase. To make such a con- 
nection, it is useful to compute the anomalous density 
matrix from the unprojected wave function defined by 
the pairing matrix that corresponds to our lowest energy 
projected state. It is given by 



Pij 



Bij 
B* ; An 



where 
An 



(*o|*o> 



Bij - 



(*ol*o) 



(3) 



(4) 



and transforms under an SU(2) gauge transformation, 
—> G • *, like pij — > G; • ■ Gj exactly like the mean 

fields in the corresponding slave particle theory (see Ref. 

[2"0] or the SM S-I). We can therefore use this matrix to 

study the projective symmetry properties of the obtained 

unprojected state |^o)- 

To study fluxes through the lattice, we can use the 

SU(2) matrix version of a "phase" variable 



= -ipij/\pij\ 



(5) 



that is an analog of the U(l) phase e MiJ we associate with 
ordinary electricity and magnetism on a lattice. Follow- 
ing Ref. [50J the analog of flux through any loop on the 
lattice is therefore the product of this matrix around the 
loop 



ijk...l 



w 



(6) 



where Ni oop is the number of bonds ij that form the 
loop. Unfortunately, this product is not gauge invariant, 
but for every loop of the lattice, we can define an angle 
8 associated with its flux through the trace of this ma- 
trix. Unlike a U(l) flux, here 6 — 2ir introduces a phase 
change of —1. The natural first loops to consider when 
characterizing the state |f2) are the nearest neighbor bow 
ties (product of two neighboring triangles) and hexagons 
(the trace of $ijk...i for an odd site loop vanishes by time 
reversal symmetry). These loops allow us to determine 
which of the four U(l) spin liquid states is closest to our 
best optimized state. The results are 



'hex, 



= (1.994 ±0.003)vr, 



'bow , 



= (0.010 ± .007)vr (7) 



where the error estimate is the standard deviation and 
(. . .) denotes the average value of the flux over differ- 
ent hexagons or bow ties respectively. These results are 
equivalent to nearly tt flux through the hexagon and 
flux through the bow tie or triangle in the traditional 
U(l) description of flux. They show distinctly that the 



optimized state is very close to the U(l) Dirac state as 
expected from energetic considerations but that because 
of the symmetry breaking shown in Fig. [T] it is an insta- 
bility. 

Having established that our newly discovered state is 
very close to the Dirac state, we now turn to its symme- 
try breaking properties. To this end, we seek the space 
group representations that give the largest contribution 
to the pattern in Fig. [ljb) which are periodic with at 
most a quadrupled unit cell. Such representations were 
constructed in Ref. [3] and labeled as At, A 2 , B\, B 2 , 
Ei, E 2 for those related to the point group alone and 
Fi, F 2 =F 1 (g) A 2 , F 3 = Fi® B u F 4 = Fx ® B 2 for those 
allowed by a doubling/quadrupling of the unit cell. The 
focus of Ref. [4] was on the F\ representation for the 
"Hastings valence bond crystal" states associated with 
the generation of mass of the Dirac fermions. However, 
the bond amplitudes plotted in Fig. [ljb) are not of this 
representation. Instead, they are dominated by the F 2 
and E 2 representations whose patterns are shown in Fig. 
[3] (a) and (b) . Remarkably, the symmetry of the F 2 pat- 
tern alone is the same as the symmetry of the F 2 and 
E 2 patterns. The E 2 pattern alone, however, has higher 
symmetry Hence, the symmetry breaking observed here 
arises uniquely from a desire to form the F 2 pattern. 

The only time reversal symmetric alternative to 
the Hastings states, among instabilities of the Dirac 
fixed point identified in Ref. [4] is the spin sin- 
glet/nodal triplet "w-monopole" that is created by 
a complex operator Wi, i — x, y, z]$\. In SM 
S-IV we show, following the transformation prop- 
erties determined in Ref. [4 , that the six dimen- 
sional vector (Re w x ,Kew y , Re w z , Im w x , Im w y , Im w z ) T 
transforms under the two three-dimensional representa- 
tions Fi and F 2 . This remarkable coincidence allows us 
to conjecture that the w-monopole is responsible for the 
instability of the Dirac state observed in Fig. [ljb) . 

To further understand the role of the symmetry break- 
ing, we measure the asymmetry of Fig. [I] during part of 
an optimization run and correlate this with the observed 
energy. We compare global asymmetry (removing the E 2 
representation to remove the effects of twisted boundary 
conditions) defined as 

O = £<*o|& ■ Sj\* Q ) - El^Si ■ £.,1*0) (8) 

as well as the assymetry of the F 2 component defined as 
Of 2 =]>> 1 lJ (vI/ |5 l -5 J |*o) (9) 

(ij) 



where F 2 \j = 1 



2Nb/3 on the solid thick bonds, F 2i j = 
— 1/ W2Nb/3 on the dashed thick bonds and zero oth- 
erwise with Nb the number of nearest neighbor bonds. 
These are shown in Fig. [3] (c) and (d). Interestingly, 
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(c) (d) 

FIG. 3. Spatial symmetry analysis of our state, (a) and (b) 
are the kagome space group representations associated with 
the pattern of Fig. [ljb) . (c) and (d) Correlation between the 
energy of different states and the amount of assymetry ((c): 
as defined by eqn. [8] (d): as defined by eqn. (9|. Red circles 
correspond to states generated during part of an optimization, 
the green star to the Dirac spin liquid and the blue triangle 
to our most optimized state. For energies below the Dirac 
spin liquid, improved energy is correlated with decreased total 
assymetry but increased strength of the Fs pattern. 

although we see a clear symmetry breaking, lower en- 
ergy states have lower asymmetry, saturating at a value 
(triangle) that is small but still above the Wen Dirac 
state (star). We find exactly the opposite behavior for 
the overlap with the F2 component: lower energy states 
have greater overlap. This implies the symmetry break- 
ing associated with the F% representation is an important 
feature of our state throughout much of the minimization 
process. 

In addition to breaking the spatial symmetries of the 
lattice, the physical system also breaks the U(l) sym- 
metry of the Dirac state down to a Z 2 symmetry. To 
study this symmetry breaking we produced a series of 
runs without any pairing in the wave functions. We ini- 
tialize twelve runs by taking the Dirac pairing function 
and multiplying each element by (1 + r) where r is sam- 
pled in the interval [0,X]. Using X = {0.1,0.2,0.33} 
we ran 12 simulations whose variation in initial energy 
was significant starting as high as E — —0.375 per site. 
All simulations converged to E = —0.4295, the energy 
of the Dirac state. From these results we conclude that 
the Dirac state lies at the bottom of a deep and wide 
region in orbital space and that the pairing of spinons is 
necessary to produce the state shown in Fig. [ljb). 

Based on our results, it seems very likely that the Dirac 
fixed point is unstable to the formation of a stripe-like 



spin-liquid crystal phase. This conclusion rests on two as- 
sumptions: that fluctuations beyond those captured by 
the projected wave function do not restore the symme- 
try and that this symmetry breaking is not a finite size 
effect. We have some indirect evidence supporting the 
latter assumption. DMRG and exact diagonalization re- 
sults indicate that a 4x4 unit cell cluster should be large 
enough to capture the qualitative physics of the system. 
In addition, we have looked at up to 8x8 unit cell clusters 
and can still find states with lower energy than the Dirac 
state. 

One way to directly address these assumptions would 
be to perform a PSG analysis on the relevant lower sym- 
metry subgroup of the kagome lattice and use it to search 
for the state that projects to our state (we cannot do 
this directly because projection is a many-to-one map- 
ping and cannot be inverted) . Such an analysis would al- 
low a determination of the finite size scaling of the sym- 
metry breaking effects and provide a starting point for 
studying fluctuations about this phase in a low energy 
effective theory. This would also establish more directly 
the question of whether there is an energy gap (however, 
we expect such a gap because the wave vector of the sym- 
metry breaking pattern connects the Dirac nodes of the 
Brillouin zone). 

Given the lack of evidence in large-scale DMRG calcu- 
lations for our state, it is unlikely to represent the true 
ground state. However, since it lies very close by and 
involves a minimal loss of crystal symmetries, it seems 
likely to be a leading instability. In particular, small per- 
turbations to the Hamiltonian could stabilize it suggest- 
ing it could be realized in nature. One promising class of 
materials are the Zn-Paratacamite family parameterized 
by Zn doping concentration x with x < 1/3. Unlike the 
structurally perfect kagome lattice of the x = 1 Hcrbert- 
smithite member of the family, compounds with x < 1/3, 
including clinoatacamite at x — 0, break crystal symme- 
tries and have distorted kagome layers [17] with precisely 
the distortion expected from the symmetry breaking of 
our state. Our results therefore motivate the study of 
single crystals of these materials and suggests that either 
the mysterious intermediate phase below 7K < T < 20K 
or the high temperature phase T > 20K could still have 
spinons as low energy excitations that are delocalized 
along the "rails" of the distorted lattice. 

The most remarkable implication of our results is 
its suggestion that spin liquid crystal phases may be 
a common phenomena. Since any dimer state is the 
exact ground state at the mean field level[5T], projec- 
tion must introduce quantum fluctuations that melt such 
crystalline phases, take the system through a succession 
of more symmetric phases until, in our case, it nearly 
reaches an isotropic phase. Such a picture has several 
implications for the DMRG calculations on the kagome 
lattice. It is known 5J that small perturbations to the 
Hamiltonian in DMRG (boundary conditions, pinning 
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fields, etc.) can enhance different states and lead to 
symmetry breaking. Exploring the class of perturba- 
tions that stabilize the symmetry breaking we observe 
here would help make a deeper connection between ana- 
lytic Schwinger-fermion theory and DMRG. More inter- 
estingly, the observed DMRG state might be understood 
as a further instability of our state to a nematic spin 
liquid crystal, a state found recently on a triangular lat- 
tice model with ring exchange [22] . A nematic spin liq- 
uid phase would be indistinguishable from a spin liquid 
phase on long cylinders that explicitly break rotational 
symmetry. Of course, it is also possible that such a puta- 
tive nematic state melts into an isotropic Zi spin liquid. 
More generically, our results suggest that spin liquid crys- 
tal states are common in frustrated antifcrromagnets and 
may be found using methods similar to ours as competi- 
tive ground states in many other systems. 
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